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ABSTRACT 

Context. Finite radius accretion disks are a strong candidate for launching astrophysical jets from their inner parts and disk- winds are 
considered as the basic component of such magnetically collimated outflows. Numerical simulations are usually employed to answer 
several open questions regarding the origin, stability and propagation of jets. The inherent uncertainties, however, of the various 
numerical codes, applied boundary conditions, grid resolution, etc., call for a parallel use of analytical methods as well, whenever 
they are available, as a tool to interpret and understand the outcome of the simulations. The only available analytical MHD solutions 
for describing disk-driven jets are those characterized by the symmetry of radial self- similarity. Those exact MHD solutions are used 
to guide the present numerical study of disk- winds. 

Aims. Radially self-similar MHD models, in general, have two geometrical shortcomings, a singularity at the jet axis and the non- 
existence of an intrinsic radial scale, i.e. the jets formally extend to radial infinity. Hence, numerical simulations are necessary to 
extend the analytical solutions towards the axis and impose a physical boundary at finite radial distance. 

Methods. We focus here on studying the effects of imposing an outer radius of the underlying accreting disk (and thus also of the 
outflow) on the topology, structure and variability of a radially self- similar analytical MHD solution. The initial condition consists of 
a hybrid of an unchanged and a scaled-down analytical solution, one for the jet and the other for its environment. 
Results. In all studied cases, we find at the end steady two-component solutions. The boundary between both solutions is always 
shifted towards the solution with reduced quantities. Especially, the reduced thermal and magnetic pressures change the perpendicular 
force balance at the "surface" of the flow. In the models where the scaled-down analytical solution is outside the unchanged one, the 
inside solution converges to a solution with different parameters. In the models where the scaled-down analytical solution is inside the 
unchanged one, the whole two-component solution changes dramatically to support the flow from collapsing totally to the symmetry 
axis. 

Conclusions. It is thus concluded that truncated exact MHD disk- wind solutions which may describe observed jets associated with 
finite radius accretion disks, are topologically stable. 

Key words. MHD — methods: numerical — ISM: jets and outflows — Stars: pre-main sequence, formation 

1, Introduction that jets are fed by the material of an accretion disk surround- 

ing the central object. The gravitational energy of the inf ailing 
Astrophysical jets are ubiquitous, occurring m a variety of ob- materM in me disk is converted int0 kinetic energy of the out _ 

jects on very different size and mass scales. They can be pro- flowing matter Radiative i aun ching can be excluded due to the 

duced by pre-mam sequence stars m Young Stellar Objects, by involved small rad iation fields which have usually luminosities 

post-AGB stars m pre-planetary (PPNs) and planetary nebulae of only a few percent of the necessary kinetic luminosities (e.g. 

(PNs), by white dwarfs (WDs) m supersoft X-ray sources and T^ a i m in the case of YSOs). Furthermore, the kinetic lumi- 

symbiotic stars (SySs), by neutron stars m X-ray Binaries, by nosity of the outflow seems t0 be a large fraction of the rate at 

stellar black holes m Black Hole X-ray Binaries and by su- which energy is re i ea sed by accretion. With such a high ejection 

permassive black holes in the case of Active Galactic Nuclei. efficiency it is natural t0 assume that jets are driven magne tically 

However, the formation mechanism of jets remains very poorly from an accretion disk; the magne tic model of a disk-wind seems 

understood. t0 ex pi am simultaneously acceleration, collimatio n as well as the 

In all such cases, jets and disks seem to be inter-related, observed high jet speeds dKonigl & Pudritzll200Ql) . 
Disks provide the jets with the ejected plasma and magnetic 

fields and jets are possibly the most effici ent means to re move In their seminal analytical work iBlandford & Pavnd (Il982l) 
excess angular momentum in the disk fe.g. lFerreirall2007l) mak- studied the magneto-centrifugal acceleration along magnetic 
ing accretion possible in the first place. Observa tionally, such a field lines threading an accretion disk. They were the first 
correlation is now already well established (e.g. lHartigan et alJ to show the braking of matter in the azimuthal direction in- 
[1995 in the case of YSOs). Hence, our current understanding is side the disk and the outflow acceleration above the disk 
surface guided by the poloidal magnetic field components. 
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iBlandford & Paynel (Il982l) along the guidelines of radially self- 
similar so lutions of the full magnetohydrodynamics (MHD) 
equat ions (IVlahakis & Tsinganoslll998h . lOgilvie & Livid (11998 . 
l200lh solved for the local vertical structure of a thin disk 
threaded by a poloidal magnetic field of dipolar symmetry. They 
showed that jet launching occurs only for a limited range of field 
strengths and a li mited range of inclination a ngles. 

The model of IBlandford & Paynel (Il982h . however, has basi- 
cally two serious limitations. First, the outflow speed at large dis- 
tances does not cross the corresponding limiting characteristic, 
with the result that the terminal wind solution i s not ca usally dis- 
connected from the disk. Later, Vlaha kiset al.l (I2QQQL V00 here- 
after) showed that a terminal wind solution can be constructed 
which is causally disconnected from the disk and hence any per- 
turbation downstream of the superfast transition cannot affect 
the whole structure of the steady outflow. The second limita- 
tion of self- similar models in general is their geometrical na- 
ture. Singularities exists at the jet axis in radial self-similar mod- 
els. Numerical simulations are necessary to extend the analyti- 
cal solutions, as it has been already done in [Gratia et al.l (120061 
GVT06 hereafter) and lMatsakos et alJ <|2QQ8 M08 hereafter). 

Another geometrical limitation of radially self- similar mod- 
els is the non-existence of an intrinsic scale. The jets formally 
extend to radial infinity. The aim of this paper is to investigate 
numerically, how imposing an outer radius of the jet, i.e. cutting 
off the analytical solution at arbitrary radii, affects the topology, 
structure and stability of a particular radial self- similar analytical 
solution and hence its ability to explain observations. 

Several other numerical studies exist, which have focused 
on the magnetic launching of disk- winds. In most models a 
polytropic equilibri um accretion disk was regarded as a bound 
ary c o ndition (e.g 



V B = 0, 



(5) 



19991 120031: lOuved 



Ustvusova et al. 
detail l2003p 



skv et al. 

Anderson et aDl2005l 120061: iPudritz etaLll2006h . The magnetic 



19951: iKrasnopolsky et al.l 



Nakamura & Meie 



feedback on the disk structure is therefore not calculated self- 
consistently. Only in recent years were the first simulations 
including the accretion disk self-con sistently in the calcula- 
tions of jet formation presented (e.g. ICasse & Keppensl [2002, 
12004 iKato et al.l 12001 IZanni et al.l 120071) . Nu merical simula- 
tions of stellar winds have also been done by e.g. lMatt & Pudritzl 
However, in none of these studies the disk has been 
truncated to study the effects of an intrinsic scale. 

This paper is organized as follows: in Sec. [2 the analytical 
self- similar model underlying our numerical study is reviewed. 
The numerical setup is described in Sec. [3] The study starts with 
several test simulations which are described in detail and whose 
results are presented in Sec. ffl In Sec. \5\ is described the pa- 
rameter study of steady solutions and its results. We close with 
a summary of the results in the last section. 



2. Analytical self-similar model 

The ideal time-dependent MHD equations which are solved nu- 
merically are: 



dp 

dt 



dv 

Hi 



P 

dp 
dt 



f x (V x B) + - V p 
P 

+ v-V/? + I>Vv = A 
dB 

— - V x (v x B 

dt 



0, 


(1) 


-vo, 


(2) 


A, 


(3) 


0, 


(4) 



where p, p, v, B denote the density, pressure, velocity and mag- 
netic field over V47T, respectively. O = -Q Mir is the gravita- 
tional potential of the central object with mass At, A represents 
the volumetric energy gain/loss terms (A = [T - l]pg, with 
Q the energy source terms per unit mass), and Y is the ratio of 
the specific heats. The spherical radius is denoted by r and the 
cylindrical radius by R. 

By assuming steady-state and axisymm etry, several con - 
served quantities exist along the field lines (lTsinganoslll982l) . 
By introducing the magnetic flux function A = (1/2 tt) f B p • dS 
to label the iso-surfaces that enclose constant poloidal magnetic 
flux, then these integrals take the following simple form: 



¥a(A) = 



pv P 



^(A) = - v 

R \ p 



L(A)=R\v,-^-), 



(6) 
(7) 
(8) 



where *Fa is the mass-to-magnetic-flux ratio, Q the field angular 
v elocit y, and L the total specific angular momentum. The ratio 
^L/£l defines the Alfvenic lever arm Ra at the point where the 
flow speed is equal to the poloidal Alfvenic one. At the refer- 
ence field line (see below) we set Ra\o=i = ft*. In the adiabatic- 
isentropic case where A = 0, there exist two more integrals, the 
ratio of total energy flux to mass flux E and the specific entropy 
<2, which are given by: 



E(A) 
Q(A) 



v 2 r 

2 + fT 
P_ 



1 P 



+ &-Q.R 



(9) 
(10) 



We use the steady, radially self- similar solution which is de- 
scribed in V00 and crosses all three critical surfaces. We note 
that a polytropic relation between the density and the pressure is 
assumed, i.e. P = Q(A) p r , with y being the effective polytropic 
index. Equivalently, the source term in Eq. © has the special 
form 



A = (T-r)/>(V-v), 
transforming the energy Eq. ([3]) to 



dp 

+v • V p + y pV v 

dt 



0. 



(11) 



(12) 



The latter can be interpreted as describing the adiabatic evo- 
lution of a gas with ratio of specific heats 7, whose entropy P/p y 
is conserved along each streamline. 

The solution is provided by the values of the key functions 
M(6), G(6) and iff (6), which are the Alfvenic Mach number, the 
cylindrical distance at any point on a field line in units of the 
corresponding Alfvenic lever arm and the angle between a par- 
ticular poloidal field line and the cylindrical radial direction, re- 
spectively. 6 is the colatitude in a spherical coordinate system. 
Then, the poloidal field lines can be labeled by 



B*Rl 



W2 



with 



R 2 



rI g 2 



(13) 
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The physical variables can be constructed as follows (V00) 
1 

1 



p* a 



,x-3/2 



p* a' 



x-2 



M 2 y ' 

M 2 



v D = -v* a 



= v*A a' 



-1/4 



sin# 



-1/4 



G 2 cos(<A + 6>) 
G 2 -M 2 
G(l - M 2 ) ' 



(cosiJ/e R + smi^e z ) , 



(14) 
(15) 
(16) 
(17) 



B p = -B* a x/2 ~ l S1 ^ (cosi/re R + sinif,e z ) , (18) 

G z cos(^ + 6) 



B* ~ -B*Aa 



,x/2-l 



l-G 2 



G(\- M 2 ) ' 



(19) 



Note that x is a model parameter governing the scaling of the 
magnetic field and is related to £ = 2 (x - 3 /4) which i s a local 
measu re of the disk ejection efficiency in the model of iFerreiral 
(1 19971) . The starred quantities are related to their characteris- 
tic values at the Alfven radius R* along the reference field line 
a = 1 . Moreover, they are interconnected with the following re- 
lations: 



B* 



2 



7C = 



\QM 

v? 



(20) 



The constants A and 7C are the specific angular momentum of the 
flow in units of v* R* of the reference field line and the Keplerian 
velocity at distance R* measured on the disk in units of v*, re- 
spectively. Finally ft is the ratio of gas pressure to magnetic pres- 
sure and is proportional to the gas entropy. 

3. Numerical setup 

We solve the MH D equations with the PLUTC0 code 
dMignone et al. I [2007b . a modular Godunov-type code particu- 
larly oriented towards the treatment of astrophysical flows in the 
presence of discontinuities. For the present case, second order 
accuracy is achieved using a Runge-Kutta scheme (for temporal 
integration) and piecewise linear reconstruction (in space). All 
the computations were carried out with the simple (and compu- 
tationally efficient) Lax-Friedrichs solver. 

We define the reference length R* to be unity, while the ref- 
erence velocity is normalized by setting v* = 1 . Time is given in 

units of to = 2 n ^RllQM, i.e. one Keplerian orbit at R* = 1. 
The model parameters of this solution were chosen as x = 0.75 
and y - 1.05, while the solution parameters are given to be 
A = 11.70, p = 2.99, % = 2.00, corresponding to the solu- 
tion in V00 crossing all critical points. At the symmetry axis, 
the analytical solution is modified as described in GVT06 and 
M08. 

To study the influence of the truncation of the analytical so- 
lution, we divide our computational domain into a jet region and 
an external region, which are separated by a truncation field line 
^tmnc- For smaller values of a - i.e. smaller radii - our initial 
conditions are fully determined by the solution of V00. The way 
the external region should be initialized, however, is not as obvi- 
ous. 

Therefore the technical questions, which arise at this point 
and which are tested in several simulations described below, are 
the following: 
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- What are the consequences of different initializations of the 
external region? 

- How is the choice of canine affecting the results? 

- What are the effects of numerical resolution and domain size. 

We impose axisymmetry at the axis, i.e. variables are sym- 
metric across the boundary, normal components of vector fields 
and also angular components change sign there. Outflow condi- 
tions are set at the top z and outer radial boundaries, i.e. all gradi- 
ents across these boundaries are set to zero. At the lower bound- 
ary, we keep the quantities fixed to their analytical values, how- 
ever, making sure that the problem is not over- specified. Thus, 
the number of quantities set by the analytical solution is lowered 
by one for each critical surface already crossed upstream of the 
boundary (for details, see M08). At small radii, where the flow 
is super-fast-magnetosonic, seven of the eight MHD quantities 
are given by the analytical solution, vr is set by the require- 
ment vp || B p . At larger radii, where the flow is sub-fast and 
super- Alfvenic, only six quantities are fixed, is set by sym- 
metrizing the values inside the domain across the boundary. In 
the sub-Alfvenic and super-slow regime, also p is given by its 
values inside the domain, and in the sub-slow regime, also B R is 
given. 

4. Test simulations 

4.1. Initial conditions 

To answer the questions posed above, we have run ten models, 
together with the unchanged model of V00, hereafter labelled 
ADO (analytical disk outflow solution, as in M08). Details of 
these test models are given in Table Q] In all these models the 
internal part of the flow is initialized with the V00 solution. 
We probe the second question by comparing the models FL1 
- FL3, where the analytical solution is truncated at a tmnc = 0.4, 
#trunc =0.8 and a tmnc = 5, respectively. In these models the 
external region is initialized with constant values, the values of 
the analytical solution at the point (R = R max , Z = Z max ). The 
influence of numerical resolution and domain size is tested with 
the models CD1 - CD5, in which the external region is initial- 
ized using the analytical solution, but with p and v z damped with 
^tnmcM- The first and most important question can be addressed 
with models FL1, CD1 and ER1 - ER2. In model ER1 we use 
the analytical solution in the external domain, but v z is set to a 
smaller value. In model ER2 all quantities are damped with an 
exponential factor in the external domain. 

4.2. Results of the test simulations 

The basic evolution is similar in most models - as an example 
we show the structure of the flow for models CD1 and ER2, see 
Figs. [T]and [2]) - and can be divided into four phases: 

1. At the beginning, a shock front starting at the jet base runs 
across the jet, bending its outer surface and forming a dent 
which then travels out- and upwards (see Fig. [2 third panel, 
black lines). The front consists of two distinct shocks, pre- 
sumably slow- and fast-magnetosonic waves which transport 
downstream the effect of the boundary condition at the base, 
namely the truncation of the solution. They can be seen in all 
models, but are not present in the analytical reference model 
ADO as expected. A third feature, which is also present 
in model ADO, is the fast magnetosonic separatrix surface 
(FMSS) which shields the flow from the modification close 
to the axis (GVT06, M08). 
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Table 1. List of numerical test models 



Name 


Grid size 


Resolution 


^end 


*coll 


Description 


ADO 


[0,50] x [6,100] 


200 x 400 


50.0 




model of Vlahakis et al. (2000), Gracia et al. (2006) and Matsakos et al. (2008) 


Models to test the influence of different truncation field lines 


FL1 


[0,50] x [6,100] 


200 x 400 


8.8 


5.7 


^trunc = 0-4, constant values in external region (analytical values at LR m ax,Z max ]) 


FL2 


[0,50] x [6,100] 


200 x 400 


7.3 


6.7 


<^trunc = 0.8, constant values in external region (analytical values at LK m ax,Z max ]) 


FL3 


[0,50] x [6,100] 


200 x 400 


10.0 


8.1 


QWc - 5., constant values in external region (analytical values at LR m ax,Z max ]) 


Models to test the influence of the computational domain 


CD1 


[0,50] x [6,100] 


200 x 400 


6.3 


2.7 


tftrunc = 0.4, in external region damping of p and v z with a tmnc /a 


CD2 


[0,100] x [6,200] 


200 x 400 


8.9 


3.8 


same as CD1, but with larger domain (lower resolution, same grid size) 


CD3 


[0,100] x [6,200] 


400 x 800 


8.1 


4.1 


same as CD1, but with larger domain (same resolution) 


CD4 


[0,200] x [6,400] 


400 x 800 


9.8 


6.9 


same as CD2, but with even larger domain 


CD5 


[0,400] x [6,100] 


1600 x 40C 


1 4.7 


4.5 


same as CD1, but with eight times larger domain in R direction, same res. 


Models to test different initializations of the external region 


ER1 


[0,50] x [6,100] 


200 x 400 


50.0 


24.9 


Q^trunc = 0.4, in external region damping of v z to constant value 












(10% of the value at the truncation field line) 


ER2 


[0,50] x [6,100] 


200 x 400 


3.6 




owe - 0-4, in external region damping of all variables with exp[-(af/af trun c) 2 ] 



2. Behind the dent, a new smooth jet surface develops, some- 
times with a larger radius than the initial one. This configu- 
ration is stable for several to in all models. 

3. In the next phase, the outer medium prevails and compresses 
the jet along its full length, pinching it even close to the 
rotation axis at certain z values. Here very dense knots are 
formed. We call this phase the collapse of the solution, since 
the jet width becomes smaller than the outer disk radius (e.g. 
5.375, if a t runc = 0-4). The only model without this collapse 
is model ER2. 

4. Finally, jet material pushes back and a more or less cylindri- 
cal topology is established. The jet pulsates along the whole 
jet length in the computational domain and the time step de- 
creases by three up to nine orders of magnitudes, making it 
impossible to further advance the simulations. 

The fact that a new smooth jet surface forms right behind 
the dent is a promising result, which may be a hint for a new 
equilibrium and thus a steady solution. However, this couldn't 
be tested in these runs, since the solution collapsed. 

In Fig. [3]we plot the radii of the jets in models CD1, ER2 
and ADO at several values of z- These radii correspond to the 
truncation field line (anchored in the lower boundary where a = 
^tmnc at R = 5.375), although this is not the outermost field line 
inside the jet for all timesteps in all models. 

In model CD1, the jet radius stays constant until the collapse. 
The collapse starts at high latitudes (about z = 50) and then af- 
fects the jet further down-stream. After the collapse, the radius 
varies along the jet length and also with time. A similar structure 
is also seen in models CD2-CD4 and FL1-FL3 with similar val- 
ues of the jet radii. 

The time t co \\ at which the jet collapses is different in each 
model (Table [T]). In models CD1-CD4 testing the influence of 
numerical resolution and domain size we can clearly see that the 
size of the computational domain seems to be important for trig- 
gering the collapse. The larger the computational domain, the 
later the collapse occurs. The similarity between models CD2 
and CD3 shows that the resulting timescale does not depend on 
the numerical resolution, or that the coarsest numerical resolu- 
tion we chose seems to be enough to properly describe the prob- 
lem. 

It is very likely that the collapse is a numerical artefact trig- 
gered by the boundary conditions at the outer radial boundary. 
Since we used "outflow" conditions, all gradients are set to zero 




Fig. 3. Evolution of the jet radius (radius of the truncation field 
line) for models CD1, ER2 and also the unchanged analytical so- 
lution ADO, at two representative z values above the equatorial 
plane (z = 10 and z = 25 in the models ER2 and ADO; z = 10 
and z = 50 in the model CD1). The models ADO and ER2 are 
stationary while the model CD1 collapses. 



at the boundary. This choice affects the radial component of the 
Lorentz force 



l d J± 

~2 dR 



R 



1 dB 1 

2 ~dR 



z dZ 



(21) 



i.e. it cancels the radial gradient of the magnetic pressure, but it 
leaves the pinching force (second term on the right-hand side) 
unaf fected. This may result in artificial collim ation of the flow 
(e.g. lUstvugova et al.l [T999I: IZanni etaDl2007h . Only in model 
ER2, where the toroidal magnetic field component B<p (not only 
the density and the Z velocity component) is damped to zero at 
the outer radial boundary, the collapse can be avoided. Although 
the simulation time is very short (3.6 to), we can see that the 
model is stationary up to this point, after a short period of ex- 
pansion due to higher thermal and magnetic pressure in the in- 
ner region which is also present in the science simulations (see 
below). 

Other effects caused by the very small values of p and vz in 
the external region as in model CD5 may be present. The fact, 
whether the outer radial boundary is causally connected with the 
jet from the very beginning or not, seems to be irrelevant for the 
presence or absence of the collapse. 
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Fig. 1. Structure of the flow (logarithmic density plots) for model CD1 at timesteps t = 0, 2, 3, 6 to, respectively. This is an example 
for the collapse occurring in all test runs except for model ER2. The truncation field line (the magnetic field line with a = a tmnc - 0.4 
which is anchored in the lower boundary at R = 5.375) is also plotted (white line). 



-2.3 




20 30 
R 



20 30 
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Fig. 2. Structure of the flow (logarithmic density plots) for model ER2 at timesteps t = 0, 1, 2, 3 to, respectively. This is the only 
model without a collapse. The white line is the truncation field line (the magnetic field line with a = a tmnc - 0.4 which is anchored 
in the lower boundary at R = 5.375). The positions of the slow- and fast-magnetosonic waves emitted by the truncation field line 
and of the fast magnetosonic separatrix surface (FMSS) are marked by black lines. 



In conclusion, it seems to be imperative that either the 
toroidal magnetic field component B<p should be very small at 
the outer radial boundary or that all quantities should be modi- 
fied self- consistently in order to maintain an equilibrium in the 
external region. 



In all test models, a practical problem occurs: the time step 
decreases by three up to nine orders of magnitudes, making it 
impossible to further advance the simulations. This happens in 
the collapsing models as well as in model ER2. In the latter, all 
quantities are exponentially damped to very small values close to 
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the outer radial boundary, the large gradient in the last grid cell 
meets the zero gradient boundary condition in the ghost cells, 
which triggers leaking of matter into the domain and creating 
large velocities of all waves close to the boundary. This results 
in very small timesteps. In order to continue our studies and start 
a parameter study we are thus forced to "truncate the truncation" 
as described in the next section. 



5. Parameter study of steady solutions 

5.1. Initial conditions 

In the science models SCI a - SC5, we use the following ap- 
proach: we modify again all quantities in one of the two compo- 
nents of the outflow, and second, we initialize the external region 
with another analytical solution with slightly changed parame- 
ters. One can show that if [p (R, Z) ,p(R,Z),v (R, Z) , B (R, Z)] 
is a solution of the MHD equations then 
\p' (R\Z') ,p' (R\Z') ,v' (R' ,Z') ,B' (R' ,Z')\ with 



R' = MR, 
B' =A 2 B, 



Z' = AiZ 



At' 




(22) 



is also a solution of the same set of equations. These transforma- 
tions have also implications for the integrals of motion 



Q'(A)= J-jQ(A), 



X ¥ A '(A) = A 2 /-iy A (A), 
L'(A) = VAi A 3 L(A) , 



G'(A) 



E'(A)=-^E(A), 



-,2-21 



r 



Q(A). 



(23) 



Since in our case, both solutions should have the same central 
object, we set A3 = 1. We also ignore the scaling of R and Z, 
i.e. we do not initialize e.g. p ' (7? ', Z '), but p ' (R, Z) in the outer 
region. Formally, this is not a solution of the MHD equations, 
since the gravity term explicitly depends on the length scale. In 
our case, however, the gravitational force is small compared to 
the other forces, thus this slight inconsistency is unimportant. 
We match both solutions by using the function 

Q = Qjet exp[-(a/aw) 2 ] + Qext (l - exp[-(a/a tmnc ) 2 ]) (24) 

for all quantities. 

Note that model ER2 can be also described by this ansatz, 
when we set Q QXt = , i.e. simply damp all quantities with the 
exponential factor, or equivalently if we set A 1 — > 00 and A 2 — > 
0. However, this truncation led to problems with the time step 
(see Sec. 01). In model SCla, we try to mimic this behavior less 
drastically by setting A\ = 10 3 and A 2 = 10" 3 (v' = 10" 3/2 v, 
B f = 10~ 3 B,p' = 10" 6 p,p' = 10" 3 p). 

In Eqs. (f22l)-(f23lh several special cases can be distinguished, 
coinciding with combinations of A\ and A 2 where quantities or 
integrals remain unchanged in both solutions: 

1. If A\ - \,R,Z,v, Q(A), L(A) and E(A) remain the same. 

2. If A 2 = 1, B and p are the same. 

3. If yfJi A 2 = 1, p and *¥a (A) remain the same. 



The additional case, where A\ - A 2 - 1, is of course identical 
with the unchanged analytical solution of V00 and therefore un- 
interesting in this study. For each choice of A\ and A 2 we can cre- 
ate two different models, depending on whether the unchanged 
solution is inside or outside the solution with primed quantities. 

In model SC2, we use the scaling parameters A\ - 100 and 
A 2 = 0.1 (case 3.), Lev' = 0.1 v, B ' = 0AB,p' = 0.01 p, 
p' - p, in the exterior region. In model SC3 we use the same 
scaling parameters, but apply the primed solution in the interior. 
In model SC4, our combination is A\ = 1 and A 2 = 0.1 (case 
1.) in the exterior and in model SC5 in the interior, leading to 
v' = v, B ' = 0AB,p' = 0.01 p,p' = O.Olp. 

5.2. Results of the science simulations 
5.2.1. Structure of the flow 

As expected a priori, the flows behave qualitatively very differ- 
ently, depending on whether the scaled-down solution is inside 
or outside the original one. 

In the cases where the quantities are scaled down in the ex- 
terior region, the opening angle of the flow increases, see Figs.@] 
and[5]). The jet surface is pushed outwards by the higher thermal 
and magnetic pressure in the inner region, which also dilutes 
while expanding. A new equilibrium is established within sev- 
eral to and is stable for at least 250 to. The final opening angles 
seen in the density plots seem to be too large for calling the flow 
a collimated jet. Using the truncation field line plotted in the fig- 
ures, we find angles of about 40°-50° in models SCla - SC2 
and model SC4. In paper II of this series, we assess the "appar- 
ent" opening angle by calculating emission maps and find that 
the emitting region is very collimated in all models (the opening 
angles are of the order of 10°-20°). 

If all quantities are reduced in the internal region (models 
SC3 and SC5), then the opening angle decreases. Again a new 
equilibrium is established which is stable for at least 250 to, see 
Fig. [5j The final opening angles in these runs are very small, 
around 5°. New features are several shocks, the innermost of 
which is also collimated following magnetic field lines. 

The expansion of the flow in models SCla-SCle, SC2 and 
SC4 (Fig.©, and the collimation in models SC3 and SC5 (Fig.© 
is clearly visible in the evolution of the jet radius. While the jet 
radius is very similar in both models showing the collimation, a 
clear trend can be identified in the first case: the largest expan- 
sion occurs in model SC4, followed by the models SCla, SClb, 
SC2, SClc, SCld and SCle. Surprisingly model SCla, which 
is the model with the smallest values of thermal pressure and 
magnetic field in the outer region, is not the one with the highest 
expansion. Although models SC2 and SC4 have the same scal- 
ings for thermal pressure and magnetic field, they show different 
degrees of expansion. 

In models SC3 and SC5, again two shocks are present which 
move outwards. When the inner region shrinks due to its lower 
thermal and magnetic pressure compared to the exterior, two 
waves propagate towards larger radii emitted from the truncation 
field line. One wave travels faster than the other, and since both 
waves compress the medium, they seem to be again a slow- and 
a fast-magnetosonic wave. The fast wave has reached the outer 
radial boundary at 50 to, only its part at lower latitudes (z < 20) 
is still visible inside the domain. The slow wave stops when its 
lower part (at z = 6) reaches R - 22, which is exactly the lo- 
cation where the slow-magnetosonic critical surface crosses the 
lower boundary. Since the slow-magnetosonic wave cannot cross 
the corresponding critical surface, it is then attached to this point 
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Name Grid size Resolution t end Description 

"SCla [0,50] x [6,100] 200 x 400 25O0 a tmnc = 0.4, external analytical solution A } = 10\ A 2 = 10 -3 

SClb [0,50] x [6,100] 200 x 400 250.0 Qf trunc = 0.2, external analytical solution A x = 10 3 , A 2 = 1(T 3 

SClc [0,50] x [6,100] 200 x 400 250.0 Qf trunc = 0.1, external analytical solution A x = 10 3 , A 2 = 10" 3 

SCld [0,50] x [6,100] 200 x 400 250.0 Qf trunc = 0.01, external analytical solution A x = 10 3 , A 2 = 10" 3 

SCle [0,50] x [6,100] 200 x 400 250.0 aw = 0.001, external analytical solution A x = 10 3 , A 2 = \0~ 2 

SC2 [0,50] x [6,100] 200 x 400 250.0 aw = 0.4, external analytical solution A x = 100, A 2 = 0.1 

SC3 [0,50] x [6,100] 200 x 400 250.0 same as model SC2, but solutions are swapped 

SC4 [0,50] x [6,100] 200 x 400 250.0 ow = 0.4, external analytical solution A x = 1, A 2 = 0.1 

SC5 [0,50] x [6,100] 200 x 400 250.0 same as model SC4, but solutions are swapped 



and develops a standing and steepening shock inside the domain. 
A third shock, which is also present in the model ADO, devel- 
ops at the fast magnetosonic separatrix surface (FMSS). This 
weak shock acts as a "wall" protecting the sub-fast flow from 
the imposed modifications close to the axis (GVT06, M08). Note 
that the poloidal field lines develop a bend at the position of all 
shocks. 

In models SCla and SClb, a small dip is present after about 
70 to and 150 to, respectively. Before this dip, the radii are almost 
constant, as well as in the other models (Fig.[6j top, dashed line). 

5.2.2. Are there boundary effects in the science runs? 

We have also run models SCla', SC3' and SC4' with a larger 
domain ([0, 100] x [6,200]) to check whether the science runs pre- 
sented above are affected by the boundaries. 

The results concerning the jet radii measured with the trunca- 
tion field line are only changed by a few percent in models SCla' 
(with respect to the constant part in model SCla) and SC4' , with 
less effects in model SC4' than in model SCla', see Fig. [6) The 
small dip seen in model SCla is smoothed out, thus it was in- 
deed created by boundary effects as suspected. In model SC3', 
larger changes exist at higher z values (z > 75) after about 50 
to. While in model SC3, the radii decrease, they stay constant in 
model SC3'. 

Therefore, throughout the remainder of the paper, the results 
of all runs at 50 to will be considered as "real" final stationary 
states of the solutions. 



5.2.3. Force balance along the jet boundary 

The interplay between the R components of the pressure gradient 
and the Lorentz force (which, as expected, always dominates in 
those radially self-similar disk-wind models) is responsible for 
collimating the flow or triggering its expansion in all models (see 
Figs. 0- [8]). In the expanding models SCla-SCle and SC2, the 
Lorentz force and the pressure gradient are directed outwards, 
while they are directed inwards in the models SC3 and SC5. The 
first is also true for model SC4, the pressure gradient turns in- 
wards between z = 7 and z - 14. Furthermore, in all models 
except for SC4, the pressure gradient is only important at lower 
z values between 6-8. 



5.2.4. Quantities at the outer axial boundary and the 
properties of the jet 

In Figs. I9HT31 all eight magnetohydrodynamic quantities are 
plotted along the outer axial boundary at z - 100 for our models 
SCla - SC5. In each panel, we show the initial profile of our 
model (solid line), the profile of the unchanged analytical solu- 
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Fig. 6. Comparison of the jet radius evolution in models SCla 
and SCla', models SC3 and SC3' and models SC4 and SC4' 
(black dashed: unprimed models, red solid: primed models) 



tion (dash-dotted line) and the final profile at t = 50 to (dashed 
line). One can clearly see how the initial profile follows the an- 
alytical solution in the interior region and is damped in the ex- 
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model SCla model SClb model SClc model SCld model SCle 




10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 

R R R R R 



Fig. 4. Structure of the flow (logarithmic density plots) for models SCla - SCle (from left to right) at timesteps t = (top), t = 25 
(middle) and t = 50 to (bottom), respectively. Also plotted is the magnetic field line anchored in the lower boundary where ol — canine 
(white line)) . At t - 50 to, we also plot the two field lines with half and twice the radius of that of the truncation field line which are 
used in Sec. [5^31 Fig. [HI 
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model SC2 model SC3 model SC4 model SC5 
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Fig. 5. Structure of the flow (logarithmic density plots) for models SC2 - SC5 (from left to right) at timesteps t = (top), t = 25 
(middle) and t = 50 to (bottom), respectively. Also plotted is the magnetic field line anchored in the lower boundary where a = a tmnc 
(white line). At t - 50 to, we also plot the two field lines with half and twice the radius of that of the truncation field line which are 
used in Sec. ET231 Fig. [141 
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Fig. 7. R (top) and z (bottom) components of forces along the 
jet boundary at the final simulated timestep (at 50 ^o) in models 
SCla. The components are measured at the truncation field line 
plotted in previous figures and multiplied by R^ with Rr the local 
cylindrical radius of the field line 



terior in models SCla-SC2 and SC4 and vice versa in models 
SC3 and SC5. 

In models SCla-SCle, SC2 and SC4, at large radii, most 
profiles have a similar shape as the unchanged analytical solu- 
tion, we started with. The original break is mostly leveled out in 
all of the profiles. However, at small radii (R < 7) the final so- 
lution deviates drastically. This is expected, since the analytical 
solution is modified close to the axis to avoid the intrinsic sin- 
gularities. In density, pressure, axial (z) and toroidal (0) velocity 
components, as well as in axial magnetic field, a peak is present 
where all quantities are higher than at larger radii (R > 7). 
The radial (R) velocity and magnetic field components and the 
toroidal magnetic field on the other hand show strong deviations 
from the initial profiles with strong declines close to the axis 
in the former two quantities and a steep rise and a minimum at 
about R = 4 in the latter. The relative height of the central peak, 
the height of the remaining break and the depth of the minimum 
in B(p are different across the five models. Although the models 
SC2 and SC4 are special cases, where the density and the veloc- 
ity, respectively, is unchanged in the inner and outer region, the 
final profiles in all models including these two are similar. 

In models SC3 and SC5, the initial profiles have to be 
changed dramatically to stop the decrease of the opening angle 
and the collapse of the jet. The profiles are much more com- 
plicated and less smooth than in the other five cases described 
above. Again some final profiles show a similar shape as the 
unchanged analytical solution at large radii. However, superim- 
posed on all of them is a shock structure with jumps in density, 
pressure, radial magnetic field and toroidal velocity (at about 



Fig. 8. Same as in Fig. [71 but for model SC3. 



R = 4, R = 13, R = 25, and R = 37 in model SC3). The radial 
magnetic field changes sign at the innermost shock where the 
density and pressure rise steeply and a minimum in axial veloc- 
ity is present in model SC3. In model SC5, the radial magnetic 
field also changes sign at the innermost shock which is, however, 
now characterized by a steep decline in density and a steep rise 
in axial velocity. 

5.2.5. The integrals of motion 

In Fig. [141 we plot the integrals of motion given by Eqs. © - ©, 
along the truncation field line at the interface of both regions, 
along an inner field line which is anchored at half of the radius 
of the truncation field line, and along an outer field line anchored 
at twice the radius. 

When we compare the integrals along the inner field line in 
the models SCla-SCle, SC2 and SC4 with those in the ana- 
lytical model ADO, we can see that the flow in the inner region 
is very similar to the analytical solution we started with. An ex- 
ception is model SCle, where the inner field line is already very 
close to the symmetry axis, i.e. our modifications there affect the 
integrals of motion. All integrals of motion converge smoothly 
to an asymptotic value. The deviation from this value are within 
6 % for z values above 20. 

In the remaining models SC3 and SC5, the behavior of the 
integrals of motion along the inner field line is very different 
to the other models, but similar to each other. At about z = 50, 
peaks and dips are visible in all integrals, however, with the most 
pronounced in Q. 

In the outer region, all integrals seem to converge in models 
SCla-SCle, SC2, SC3 and SC5. The turnovers in model SC2 
are a result of a turnover in the field line itself. In model SC4, 
L and have already converged, the other integrals still vary. 
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Fig. 14. Integrals of motion *¥a (A), £2(A), L(A), E(A) and 2(A), normalized to their values at their end point, where they leave the 
domain, along the truncation field line (left), along a second field line anchored at half of the radius (middle) and along a third field 
line anchored at twice of the radius (right) for models SCla, SC2 - SC5 and ADO (from top to bottom). 



In conclusion, in most of our models also the external region 
reaches a steady state. 



5.2.6. Topology of the current lines 

In Fig. [T53 we plot the poloidal currents for our nine science 
models SCla - SC5. In model ADO, as well as in the inner re- 
gion of models SCla-SC2 and SC4, which remains unchanged 
with respect to the analytical solution, a re-adjusted FMSS is 
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Fig. 9. Quantities at the outer axial boundary in model SCI a: 
density (top left), pressure (top right), velocity components in R 
(second row, left), z (second row, right) and cp (third row, left), 
magnetic field components in R (third row, right), z (bottom left) 
and (p (bottom right); given are the initial profiles (solid), the 
final profiles at 50 to (dashed) and the profiles in the unchanged 
analytical solution of V00. 
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Fig. 10. Same as in Fig.H but for the model SC2. 



visible as shock in the density plots (Figs. |4]— B cf. GVT06 and 
M08). The currents in these models are counterclockwise up- 



Fig. 11. Same as in Fig.0 but for the model SC3. 
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Fig. 12. Same as in Fig.0 but for the model SC4. 



stream and downstream of the FMSS. Only along the FMSS the 
current close, i.e. the lines are clockwise. This distribution of 
currents, and in particular the direction of the resulting J p x 
force, is consistent with the decollimation and deceleration that 
the flow experiences as it passes through the shock. Thus, one of 
the effects of the new FMSS is to bend the streamlines away from 
the z-axis avoiding the overcollimation property of the origi- 
nal analytical solution. The collimation and decollimation pro- 
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Fig. 13. Same as in Fig. but for the model SC5. 




Fig. 15. Poloidal currents RB^ = const of models SCla-SCle 
(top) and models SC2-SC5 and ADO (bottom) at the final 
timestep (at 50 to) are plotted. 



cesses that can be derived from such a plot are also discussed in 
GVT06. Towards the outer region, the morphology of the current 
lines is distorted and of completely different topology. In models 
SC4 and SCla-SCle, the outer region is separated by a second 
shock with a current sheet from the inner region. 

In models SC3 and SC5, the FMSS is also present. Again 
at the shock, a current sheet separates two region with counter- 
clockwise currents. Furthermore, closer to the rotation axis, an- 
other current sheet forms an X shape with the FMSS at (R, z) = 
(5, 30) and a third current sheet in the lower right corner of the 
domain develops. 

6. Summary 

In this paper, we have studied the effects of imposing an outer 
radius of the underlying accreting disk, and thus also of the 
outflow, on the topology, structure and time-dependence of a 
well studied radially self- similar analytical solution (V00). We 
matched an unchanged and a scaled-down analytical solution of 
V00 and found in all these cases steady two-component solu- 
tions. 

We showed that the boundary between both solutions is 
always shifted towards the solution with reduced quantities. 
Especially, the reduced thermal and magnetic pressure change 
the perpendicular force balance at the "surface" of the flow. In 
the models where the scaled-down analytical solution is outside 
the unchanged one, the inside solution converges to another an- 
alytical solution with different parameters. In the models where 
the scaled-down analytical solution is inside the unchanged one, 
the whole two-component solution changes dramatically to sup- 
port the flow from collapsing totally to the symmetry axis. 

A result of the modification of the analytical solution at the 
symmetry axis is the formation of a weak shock at the fast mag- 
netosonic separatrix surface (FMSS) which shields the flow up- 
stream of the FMSS from the imposed modifications close to the 
axis. A result of the truncation of the analytical solution along 



some poloidal field line is the formation of two compressive 
MHD shocks: a fast shock that travels quickly downstream out- 
side of the computational domain and a slow shock which as it 
travels downstream with lower speed, it is locked precisely at 
the position of the lower computational z-boundary where the 
analytical slow surface meets this boundary. 

Our truncated disk-wind solutions are stable for more than 
50 to, i.e. several orbital periods at the truncation radius. These 
solutions may be relevant to describe observed jets, since the 
jet radii are too large in the untruncated analytic disk outflow 
solution, with respect to observed jet widths. We also provide 
all quantities at the outer z boundary which can now be used as 
boundary conditions for jet propagation studies. 

In the following paper II, we calculate emission maps corre- 
sponding to such truncated disk-models and compare them with 
observations. 
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